library(data.table)
library(magrittr)
library(parallel)
library(dplyr)
library(boot)
library(INLA)
setwd("E:\\5hmc_file\\2_5hmc_yjp_bam\\ASM")
dir1="./bayes_pvalue_beta0/"
dir2="./bayes_BF/"
id=read.table("./20201118/only.1.pattern.in.DC.and.bias.in.CCHC.unitID.txt",head=T,sep="\t")
id=as.character(id$x)
group1=c("X2B_X1T","M8_M7","M6_M5","M2_M1","M48_M47","M50_M49","M28_M27","M30_M29","M26_M25","M35_M36","M18_M17","M20_M19","M22_M21","M40_M39")
i=1
fn1=paste0(dir1,group1[i],".bayes_p.txt")
file1=read.table(fn1,head=T,sep = "\t")
file1$unitID=paste(file1$chrom,file1$position,file1$ref,file1$var,sep=":")
file1=data.frame(unitID=file1$unitID,normal_reads1=file1$normal_reads1,normal_reads2=file1$normal_reads2,tumor_reads1=file1$tumor_reads1,tumor_reads2=file1$tumor_reads2)
file1=file1[file1$unitID %in% id,]
names(file1)=c("unitID",paste(group1[i],names(file1)[2:5],sep="_"))
file=file1
for(i in 2:6){
  fn1=paste0(dir1,group1[i],".bayes_p.txt")
  file1=read.table(fn1,head=T,sep = "\t")
  file1$unitID=paste(file1$chrom,file1$position,file1$ref,file1$var,sep=":")
  file1=data.frame(unitID=file1$unitID,normal_reads1=file1$normal_reads1,normal_reads2=file1$normal_reads2,tumor_reads1=file1$tumor_reads1,tumor_reads2=file1$tumor_reads2)
  file1=file1[file1$unitID %in% id,]
  names(file1)=c("unitID",paste(group1[i],names(file1)[2:5],sep="_"))
  file=merge(file,file1,by="unitID",all=T)
}
filekp=file

ref=grep(names(filekp),pattern = "reads1")
alt=grep(names(filekp),pattern = "reads2")
file=filekp[,c(ref,alt)]

ASE=function(x)
{
  if (length(which(!is.na(x)))/4>=1){
    sel=which(!is.na(x))
    rct=x[sel]
    ref=rct[1:(length(rct)/2)]
    var=rct[(length(rct)/2+1):(length(rct))]
    if (length(sel)/4>1){
      df=data.frame(y=var,Ntrials=ref+var,x1=rep(c(0,1),length(sel)/4),x2=rep(c(0:(length(sel)/4-1)),each=2))
      #header=c(rep("y",length(sel)/2),rep("Ntrials",length(sel)/2),"x1","x2")
      #colnames(df)=header
      formula = y ~ 1 + f(x2, model = "iid") + x1
      m1 <- inla(formula, data = df, family = "binomial", Ntrials = Ntrials,quantile = c(0.005, 0.025, 0.975, 0.995))
      formula = y ~ 1 + f(x2, model = "iid")
      m0 <- inla(formula, data = df, family = "binomial", Ntrials = Ntrials,quantile = c(0.005, 0.025, 0.975, 0.995))
      b10 <- exp(m1$mlik[2] - m0$mlik[2])
    }else{
      df=data.frame(y=var,Ntrials=ref+var,x1=rep(c(0,1),length(sel)/4))
      #header=c(rep("y",length(sel)/2),rep("Ntrials",length(sel)/2),"x1")
      #colnames(df)=header
      formula = y ~ 1 + x1
      m1 <- inla(formula, data = df, family = "binomial", Ntrials = Ntrials,quantile = c(0.005, 0.025, 0.975, 0.995))
      formula = y ~ 1
      m0 <- inla(formula, data = df, family = "binomial", Ntrials = Ntrials,quantile = c(0.005, 0.025, 0.975, 0.995))
      b10 <- exp(m1$mlik[2] - m0$mlik[2])
    }
    return(b10)
  }
}
Sys.time()
rt=apply(file,1,ASE)
Sys.time()
result=as.numeric(as.character(rt))
filekp$Bayes_in_DC=result
write.table(filekp,"20201118/16656_DC6v6_bf.txt",quote=F,row.names = F,sep="\t")
